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We have obtained real i Stic fluid detonation-wave profiles. 1 

This work is summarized i n Pages 2-18. 

I 

We have studied the structure of uniaxially and hydro- 
statically-compressed solids and the transfer of energy among * 

translational and internal molecular modes. This work is described 
i n Pages 19-26. 

We have developed novel computational methods for simulating 
nonequilibrium processes using Gauss' Principle of Least Constraint. 

This work is described in Pages 27-38. 
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SUMMARY 


(1) TECHNICAL PROBLEM 

We are investigating mathematical models of detonation in dense fluids and 
in solids. We wish to obtain a fundamental understanding of the interaction of 
pressure waves, or shockwaves, with high explosives. 

(2) GENERAL METHODOLOGY 

We are solving the steady-state dynamical equations of continuun mechan- 
ics, including viscosity, conductivity, and Chemical reactions, to simulate 
detonation in fluids capable of undergoing fast exothermic reactions. We are 
solving the coupled ordinary differential equations of motion of molecular 
dynamics, for solids, to simulate the initiation process preceding detonation. 
We are examining simple models designed to describe the results of these macro- 
scopic and microscopic approaches. 

(3) TECHNICAL RESULTS 

We have obtained detonation wave profiles, using a realistic dense-fluid 
equation of State. We have compared these profiles to the predictions of the 
simplified Zeldovich-Von Neumann-Doering model. 

We have developed methods for simulating the rapid uniaxial compression of 
solids suited to molecular dynamics simulation. 

We have formulated the interatomic distribution function in solids, for 
comparison with atomistic simulations and use in kinetic reaction models. 

(4) IMPLICATION FOR FURTHER RESEARCH 

The relatively long vibrational relaxation time^ occurring in shocked 
diatomic molecules with realistic values of viscosity and thermal conductivity 
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suggest that the nonequi1ibrlum distributlon of energy is more important than 
the interaction of ordlnary transport properties with Chemical reactions. We 
therefore intend to concentrate, in the next two years, on reaction initiation 
and intramolecular energy transfers In reacting molecules in the solid phase. 











I. INTROOUCTION 


Shock and detonation waves involve atomic-scale processes whlch are 
Imperfectly understood and hard to measure. For this reason we have embarked 
on a program, using modern computational methods, of studying the processes 
through which mechanical shockwave energy stimulates the release of Chemical 
energy within a detonation wave. 

Our work entails the correlation of two points of view: continuum 
mechanics and atomic mechanics. We seek to understand the limltations of the 
continuun approach in dealing with atomic-scale problems and to develop methods 
for describing these small-scale processes in a realistic way. 

Considerable information on the equilibrium and transport properties of 
simple Systems is now availab1 e(1,2). Computer simulations of dense fluids, 
coupled with thermodynamic perturbation theories(3) based on simple hard and 
soft sphere model's, make it possible to describe the pressure-volume- 
temperature equation of State accurately, for sufficiently simple materials. 

For some "realistic" potentials, such as the Lennard-Jones 6-12 interatomic 
potential, comprehensive Computer simulation studies have established the 
density and temperature dependence of the viscosities and the thermal conducti- 
vity as well(2). Much less is known about the contribution of intramolecular 
degrees of freedom to the transport properties and equation of State. Studies 
on molecules like methane and benzene indicate that the Computer simulation 
approach is now feasible, even for molecules containing several atoms(4). 

A numerical study, comparing the continuum and atomistic versions of a 
very strong dense-fluid shockwave made use of both the equilibrium and non- 
equilibriun constitutive properties of the 6-12 potential, and indicated that 
useful conclusions could be drawn by comparing the two approaches(5). A 
natural extension of that work is to consider the possibility of coupling 












Chemical reactions with the disslpatlve processes of viscoslty and thermal 
conductlvlty whlch appear in shockwaves. That coupled comblnation is a detona- 
tlon wave. 

The simplest prototype condensed-phase high explosive is probably nitric 
oxide, NO. This molecule forms molecular nitrogen and oxygen in an exothermic 
reaction generating approximately 100 kilobars pressure(6). We set out with 
the goal of simulating the detonation of liquid NO using both the continuum and 
atomistic approaches. The continuum calculations are complete. Using an 
equation of State which describes the Lennard-Jones potential throughout the 
fluid regions of the phase diagram, and transport coefficients from the Enskog 
theory, we have generated detonation profi les as functions of the energy 
release and activation energy of the dense-fluid detonation. These calcula¬ 
tions, including viscosities and thermal conductivity, are compared to the 
Zeldovich-Von Neumann-Doering model in Section II. 

Shortly after this work was begun, a research program was initiated at Los 
Alamos, designed to measure and to model the detonation of nitric oxide. This 
program includes experimental, quantum mechanical, and dynamical groups working 
toward a common goal. Related Los Alamos research, carried out by Brad Holian, 
Galen Straub, and their co-workers, has established that the equi1ibration tlme 
between vibrational and translational temperatures in dense diatomic fluids 
(nitrogen) is relatively long(7), just as in the gas phase. These long equ11 i - 
bration times, of order 100 vibration times, suggest that the complete combus- 
tion process is impractically long for an atomistic simulation. For this 
reason we have concentrated our more recent work on solid-phase initiation and 
mode equilibrium. We are developing solid-phase calculations in which 
chemically active molecules can be studied under uniaxial and homogeneous 
















adiabatic compression. That work is descrlbed in Sections III and IV. We 
expect to devote the next two years' effort to enhancing our understanding of 
solld-phase detonation initiation from the atomistic mechanical point of view 




II. DENSE-FLUIO DETONATION WAVE STRUCTURE 

The classlcal steady one-dlmenslonal ZND (for Zeldovich, Von Neumann, and 
Doering) model of detonation-wave structure can be summarized with the help of 
Figure 1. It is assumed that the reactants can assume the equilibrium unreact- 
ed thermodynamic States shown on the lower Hugoniot curve as the result of 
shock compression. Each of these States satisfies the condition that its mass, 
momentum, and energy fluxes agree with those characterizing the initial state 
(state 0). This unreacted Hugoniot is hypothetical because Chemical reactions 
will gradually transform the reactants, at least for strong enough shockwaves, 
to hot products. 

The fully reacted products can also be described by a Hugoniot curve 
satisfying the conservation laws. This is the upper Hugoniot. Its position, 
relative to the lower Hugoniot, depends upon the heat of reaction Q. The 
Zeldovich-Von Neumann-Doering detonation wave model pictures the detonation 
process as two consecutive steps: first, the reactants are shock-compressed to 
the (Von Neumann spike) state S. Next, these hot reactants slowly undergo 
Chemical reactions, proceeding to the final (Chapman-Jouguet) state CJ. During 
the reaction phase the State point follows the Rayleigh line, satisfying the 
constraints of constant mass and momentum flux necessary in a steady wave. 
Neither the transport coefficients nor the reaction rate enter into Figure 1. 
These state functions determine the space and time scales associated with the 
compression and reaction phases of the detonation wave. 

The ZND model is an oversimplification. With the theory of 1 iquids suf- 
ficiently advanced that quantitative calculations, including nonideal effects, 
can be made, there is no real need for an inaccurate treatment. 

Molecular dynamics simulations have shown that strong dense-fluid shock¬ 
waves can be described fairly well by the Navier-Stokes equation(5). Figure 2 
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indlcates the relatively good agreement between the Navier-Stokes continuum 
calculations and the atomistic simulations. The principle disagreements are 
these: The wave is about 30% broader than the continuum predictions and the 
velocity distribution is considerably more complex than is the equilibrium 
Maxwel1-Boltzmann distribution. The alternative Mott-Smith approach(8), in 
which the velocity distribution is a weighted average of the reactant and 
product Maxwe11ian distributions, grossly overestimates the nonequilibrium 
character of the wave. 

We initially planned to carry out a similar comparison for nitric oxide, a 
simple liquid which undergoes the reaction 

2N0 = N 2 + 0 2 . (1) 

Exploratory calculations fitted to the kinetics of the reaction showed that the 
reaction zone extends for thousands of angstroms, far too long for molecular 
dynamics simulation. The liquid detonation wave prof ile that we display in 
Figure 3, together with the ZNO approximation to that profile, was calculat- 
ed(9) with an activation energy approximately half that of nitric oxide. This 
substantially reduced the spatial extent of the detonation wave. Under this 
arti f i c ial assumption there is a noticeable quantitative discrepancy between 
the simple two-step model and the solution of the coupled equations. There is 
no such dramatic difference for nitric oxide, at least for a one-dimensional 
wave. 


The instability of one-dimensional low-density detonation waves is well- 
established experimentally(10). Because this effect is at least as important 
as that of the hydrodynamic transport coefficients we see that further study of 
the one-dimensional detonation wave for nitric oxide, is not feasible, using 
molecular dynamics. 









Vlbratlonal energy transfer 1s very slow 1n dlatomlc molecules because the 
single vibrational mode is wel1-separated from the "lattice modes." The situa- 
tion is less clear cut for polyatomic molecules with many degrees of freedom. 
Can the solid-phase initiation problem be understood using equilibrium Chemical 
kinetics and equilibrium energy surfaces? Are nonequilibrium velocity di s t r i - 
butions and intramolecular energy transfer significant? The answers to these 
questions are unknown despite considerable interest and debate(ll). We are 
concentratlng our efforts on enhanced understanding of solid-phase Initiation, 
using the methods outlined in the following sections. 
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III. COMPUTER SIMULATIONS AT HIGH RATES OF STRAIN 

To exp1ore the rapid deformation associated with shockwaves in condensed 
fluids and solids we have developed special methods for solving the equations of 
motion of atoms in a rapidly-deforming f 1 u i d or solid(12). These methods were 
originally developed and applied to the simulation of viscous flows(13) and were 
later applied to solid-phase plast i c ity(14). In either case the macroscopic 
flow is described by a given strain-rate tensor Vu, where u is the hydrodynamic 
flow velocity. To incorporate the flow velocity into the equations of motion, a 
perturbation, the tensor product of Vu and Doll's tensor (the dyadic composed 
the particles' summed coordinate-momentum product qp) is added to the 
Hamiltonian function. Equations of motion result which include contributions 
proportional to the strain rate: 
q = (p/m) + q-V u ; 

p = F- V u -p; (2) 

E = - VP:V U . 

Equations (2) describe adiabatic homogeneous deformation. This scheme has 
been used to make quantitative estimates of the shear and bulk viscosities for 
dense fluids. 

Inhomogeneous compressions must be treated differently. The system is 
sheared or compressed by moving periodic images(5). Either method could be 
applied to Chemically-reacting solids. We expect to explore the homogeneous 
uniaxial compression of two Systems: The simple reaction model described in the 
next section and a "realistic" three-dimensional polyatomic model of a detonat- 
ing solid. This latter application should allow us to map out the path followed 
by energy during the equi1ibration process. 
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IV. ATOMIC DISTRIBUTION FUNCTIONS IN SOLIDS 


We are currently investigating a simple model for reacting atoms: Two 
atoms sufficlently close together react, releasing an energy Q. Equ11ibrium 
simulations (with Q=0) have been carried out to determine the density, strain, 
and temperature dependence of the "reaction rate." Figure 4 indicates the 
substantial dependence on strain biaxiality. Pastine's model(15), which 
contains an adjustable collision frequency u, can be used to fit the uniaxial 
compresslon results. 

A more fundamental approach incorporates the distribution of neighboring 
atoms. In the quasiharmonic case this can be done exactly. We have developed 
a numerical method for the partial integration of the Boltzmann factor, 
exp(-<t>/kT), over all but two particle coordinates, leading to an atom-atom 
distribution function. This function is a product of Gaussians. 

Because t calculation is somewhat intricate it is important to have a 
check available. The close-packed lattice, with nearest-neighbor interactions, 
is ideal for this purpose. Because the potential energy of a nearest-neighbor 
pair in a D-dimensional close-packed crystal Is kT/(D+l), the mean-squared 
deviation in the nearest-neighbor spacing should be 2kT/(D+l)<, where < is the 
nearest-neighbor force constant. We have verified that this result is repro- 
duced exactly by our numerical approach, and find that the convergence of 
higher-order and transverse correlations is relatively rapid and smooth as the 
number of particles in the perlodic cell is increased. In the two-dimensional 
case, for example, the ratios of the nearest-neighbor mean-squared longitudinal 
to transverse fluctuatlons for 9, 16, 25, and 36-atom crystals are 0.667, 

0.688, 0.695, and 0.696(16). 

We expect to apply this same technique to the atom-atom distribution 
function for polyatomic crystals undergoing uniax1al compresslon. This analog 
of Pastine's model’should prove useful in estimating the dependence of reaction 
rate on the macroscopic strain rate. _n_ 
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FI6URES: 


1. Diatomic reactant and product Hugoniots calculated with the Lennard-Jones 
equation of State. The well depth and collision diameter are e and o, 
respectively. The Initial volume 1s 0.8442No 3 . The initial temperature 
is 0.722e/k. 

2. Shockwave profile for a strong fluid shock carrying the triple-point argon 
to a final temperature of about one electron volt. The smooth curves are 
the Navier-Stokes solution. The points correspond to data from nonequili- 
brium molecular dynamics. 

3. Detonation profile calculated using (a) the hydrodynamic equations, includ- 
ing viscosity and heat conduction, and (b) using the ZND model. The activa 
tation energy is artifically reduced to enhance the difference between the 
two approaches. The heat of reaction is lOOe and the activation energy is 
70e. 

4. Reaction rate as a function of density for a simple binary-collision model. 
Both uniaxial and hydrostatic compressions are shown at a temperature about 
half the melting temperature. Pastine's model for the rate is shown as the 
full curve. 
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SUMMARY 


(1) TECHNICAL PROBLEM 

We are investigating atomistic mathematical models of energy transferred in 
condensed media. We wish to obtain a fundamental understanding of the inter- 
action of compressive waves with intermolecular and intramolecular energies, 
leading to the initiation of detonation. 

(2) GENERAL METHODOLOGY 

We analyze quasistatic compression by comparing equilibrium statistical- 
mechanical calculations with the results of equilibrium Newtonian molecular 
dynamics simulations. We are developing and solving noneguilibrium , equations of 
motion to describe the rapi d compression of microscopic polyatomic systems. We 
emphasize the flow of energy among the molecular degrees of freedom excited by 
the compression process. 

(3) TECHNICAL RESULTS 

We have determined the atomistic pai r distribution function for planar 
crystals undergoing both uniaxial and hydrostatic compression. We have compared 
these distributions to the predictions of the Pastine-Kamlet-Jacobs model. 

We have simulated bimolecular collisions of triatomic and hexatomic planar 
molecules (with hexanitrobenzene in mind). We analyze the normal-mode vibrations 
before and after collision and observe the transfer of energy among the transla- 
tional, rotational, and vibrational modes. 

(4) FURTHER RESEARCH IMPLICATIONS 
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I. INTRODUCTION 

Shock and detonation waves involve atomic-scale processes which are imper- 
fectly understood and hard to measure. The Systems contain too many atoms for an 
accurate quantum-mechanical treatment. Classical calculations can deal with 
Systems of hundreds or thousands of interacting atoms, over hundreds or thousands 
of vibrational periods. 

To evaluate the utility of the microscopic approach we have embarked on a 
program to develop and apply modern computational methods to the nonequi1ibrium 
dynamical processes responsible for explosive initiation. Adetailed microscopic 
view is necessitated by the violent nonequi1 ibrium nature of shock compression. 
Simulations Show that the intrinsic width of a shock front is only a few atomic 
diameters.l Thus the usual concept of equilibrium temperature, with Arrhenius 
kinetics, is unlikely to provide an understanding of the initiation of Chemical 
reactions. To achieve understanding of initiation detailed studies of rapid 
nonequi1 ibrium intermolecular and intramolecular energy transfer must be carried 
out and analyzed. 

Over the past several years there has been considerable progress in treating 
the deformation of simple monatomic fluids and solids. The techniques developed 
have been tested by intercompari$on with experimental data and with less- 
efficient older simulation techniques.^ The irreversible heating associated 
with rapid viscons and plastic deformation has been avoided by developing iso- 
thermal equations of motion.^ These equations constrain the second moment of 

O 

the velocity distribution, <(v-<v> >, to a fixed value. Analogous developments 
for treating diffusion and conduction stem from this same foundation.^ 

These new methods have not been applied to flexible polyatomic molecules. In 
order to carry out this general ization we have analyzed both quasistatic and 
rapid deformations of monatomic crystals (Section II), the rotational and c o 11 i - 
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II. DE FORMAT ION OF MONATOMIC SOLIDS 

In a series of studies we have characterized the compression of simple 
fluids and solids. We recently carried out a detailed study of equilibrium pai r 
distribution functions in deformed two-dimensional crystals.® Both Hooke's-Law 
and Lennard-Jones forces were used. The pair distributions, obtained by integra- 
ting the quasiharmoni c canonical probability density, differ from approximate 
theoretical treatments in two ways. First, contributions from al1 normal-mode 
vibrations are included. Second, reaction-rate contributions from al1 diatomic 
orientations are included. For these reasons the reaction rates follow neither 
the simple Gruneisen density dependence nor Arrhenius kinetics. Even under 
quasistatic conditions, uniaxial compression generally yields significantly 
higher collision (or reaction) rates than does the corresponding hydrostatic 
compression. 

Earlier deformation studies, under both shockwave and homogeneous condi¬ 
tions, have demonstrated the utility of incorporating the macroscopic velocity 
gradientVu directly into microscopic equations of motions. The shockwave 
studies emphasized the nonequi1 ibrium velocity distribution in the shock-front. 
Shockwave structure, bu 1 k and shear viscosities, and solid-phase high-strain-rate 
yield stengths^ have a 11 been determined in this way. Adiabatic deformation 
can be modeled by adding an extra contribution to the atomic velocities, 
q = q-7u, and a corresponding force p = -Vu»p. The two modi fi cati ons lead to 
the thermodynamic identity 

'E = -VP:7u, 

where E is internal energy and P is the pressure tensor. This adiabatic deforma¬ 
tion can be made isothermal instead by including an additional collective force 
-^p, with'J chosen to keep <p^/m> constant. 
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III. ENERGY TRANSFER IN POLYATOMIC MOLECULES 


The simplest polyatomic molecule exhibiting energy transfer is a planar 
triatomic molecule, with three vibrational degrees of freedom. In th i s case the 
two degenerate vibrational modes are coupled by Coriolis forces. Thus rotation 
of the molecule, even very slowly, leads to a geometric coupling of two vibra¬ 
tional amplitudes. 

Regular hexatomic planar molecules have nine vibrational degrees of free¬ 
dom, with four pairs of degenerate modes in addition to the symmetric breathing 
mode. At reasonable rotational velocities the coupling among the four pairs is 
negligible. Even at vanishing velocities Cordiolis coupling again causes the 
cycling back and forth between degenerate mode amplitudes. 

We have simulated bimolecular collisions of both triangular and hexagonal 
planar molecules, have analyzed their normal-mode vibrations throughout the 
collision process, and have observed energy transfer among al1 the translational, 
rotational, and vibrational degrees of freedom. Presently we are characterizing 
the energy distribution during a homogeneous compression. 











IV. THERMODYNAMICS OF RAPI D DEFORMATION 


The virial theorem relates the pressure tensor to interatomic forces. In 
the familiar monatomic form® each Newton's equation of motion is multiplied by 
the corresponding particle coordinate, summed and averaged, to give 

iF i j - VP £&rr). - (mrr) i . 
which is equivalent to 

PV = NkTI +21 r i j Fi j 

where N i s the number of atoms i n the volume pairs. For polyatomic molecules i t 
i s convenient to multiply each atom's equation of motion by the center of mass 
coordinate for its molecule. The force terms in each molecule cancel, with the 
result: 

PV = NkTI +2>ij F ij 

where N i s now the number of molecules and F.j i s the vector sum of forces of 
all atoms in molecule I due to atoms in molecule J. 

Now consider two kinds of adiabatic deformation. In "atomic" deformation 
each atom undergoes a displacement proportional to its location 

r = r-Vu 


and a corresponding force 

*p = -P'V u. 

r 

The adiabatic first-law identity, E = -VP-Vu, follows. In "molecular" deforma¬ 
tion each atom undergoes a displacment proportional to the center-of-mass 
coordinate R: 

r = R > V u 


as well as a corresponding force 

P = -P(*)‘ 7u. 

» 

The identify E r -PV:Vu again follows, provided that P is the "molecular" 
pressure. 
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SUMMARY 


(1) TECHNICAL PROBLEM 

We consider atomistic descriptions of matter, far from equilibrium, in order 
to understand the mechanisms underlying detonation i n high explosives. We focus 
on solid hexanitrobenzene. 

(2) GENERAL METHODOLOGY 

We develop equations of motion incorporating both thermal and mechanical con- 
straints. Adiabatic, isothermal, and isoenergetic processes can be simulated by 
solving these equations numerically. We study many special cases involving gases, 
liquids, and solids. 

(3) TECHNICAL RESULTS 

We have established the usefulness and the limits of "Gauss' Principle of 
Least Constraint" in classical many-body simulations. We have applied three di f- 
ferent types of feedback control to atomistic Systems in studying fluid and solid 
heat flow. We have considered the size-dependence of constitutive properties. We 
have developed a simple stress-free model for crystalline hexanitrobenzene. 

(4) IMP LI CATION FOR FURTHER RESEARCH 

The tools developed facilitate simulations far from equilibrium. The result 
i s enhanced correlative and predictive capabilities for systems undergoing nano- 
second to picosecond thermal and mechanical processes. This is a rapidly-changi ng 
field fundamental to the understanding of complex material behavior. 
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I. INTRODUCTION 


Systems far from equil1brium can only be simulated with dynamic equations of 
motion incorporating nonequilibrium temperature, velocity, and stress distribu- 
tions. During the past decade a variety of nonequi1 ibrium simulation techniques 
have been developed. During the past three years a consolidation and simplifica- 
tion has taken place: i t has been realized that simple feedback analysis can be 
applied to hydrodynamic and solid mechanics problems in a way consistent with 
macroscopic constitutive analyses. This fi e1d has been reviewed recently^*2.3 < 
The highlights are as follows: 

(1) There is an old mechanical principle, due to Gauss, "Gauss' Principle of 
Least Constraint", which is more general than Newtoni an or Lagrangian 
mechanics. Gauss' Principle makes it possible to develop equations of 
motion with nonholonomic constraints. Such constraints include, for 
instance, the restrictions that a particular process be carried out with 
a given time and space variation of temperature or energy or stress. 
Gauss' Principle has been widely tested, during the past three years, 
and found to be useful in studying steady flow far from equilibrium. 

(2) There is a new idea, due to Shuichi Nosfe(Yokahama), which makes it possi' 
ble to simulate, with molecular dynamics, systems described with indepen 
dent variables such as temperature and stress, as opposed to energy and 
volume.Nos^'s ideas are extremely interesting for theoretical 
mechanics and suggest new practical techniques, different to Gauss', for 
treating systems far from equilibrium. 

(3) Both Gauss' and Nosfe's ideas can be viewed as special cases of Control 
Theory^. This realization has increased the multiplicity of algorithms 
available, but with the advantage of a systematic description. The 












consequences of this development are now only beginning to be exp1ored. 
This is a most exciting development in mechanics, with implications for 
quantum as well as classical Systems. 

Ouring the past three years we have contributed to the development of these 
ideas while carrying out an investigation of solid hexanitrobenzene. The complex- 
ities associated with modelling this material precluded a complete study of its 
response to shockwave initiation and detonation, but the groundwork and the tools 
necessary to this simulation are now in place. 

The following sections describe the new techniques available for atomistic 
simulations, a simple model for crystalline hexanitrobenzene, and recommendations 
for further work. 
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II. NEW TECHNIQUES IN SIMULATIONS FAR FROM EQUILIBRIUM 


In 1972 Ashurst simulated the flow of momentum p and the flow of heat between 
atomistic reservoirs. He performed a momentum rescaling 

p’ = (1 + e 2 ) p + e 1 

i n which the small numbers e-| and e 2 were chosen to restore the first and 
second moments of the velocity distribution to desired values. His were the 
first systematic studies of nonequilibrium Systems using nonequilibrium molecular 

Q 

dynamics . 

In 1982 Evans and Hoover noted that the rescaling of momenta could be carried 
out i n a continuous way by adding a constraint force to the equation of motion for 
each particle: 


F c = - 2P * F 0 . 

Hoover pointed out that this same constraint force can be derived from Gauss Prin- 
ciple of Least Constraint. By 1985 i t became clear that many forms of nonequilib¬ 
rium constraint forces are possible. These were applied to a variety of Systems 

3 

far from equi11brium . 

The theoretical bases of these constraint forces advanced too. It was pos¬ 
sible to show that exact linear transport coefficients could be generated in 
steady nonequi1 ibrium simulations. Even nonlinear effects, such as the decrease 
of viscosity with strain rate and the increase or decrease of conductivity with 
heat flux, depended only weakly on the particular constraint force used. In the 
present section we summarize the two successful approaches which emerged from this 
work and point out their usefulness in studying systems far from eguilibrium. 






A. Gauss' Principle 


1 

Gauss stated that any conceivable constraint on a dynamical system should be 
implemented by minimizing the mean squared value of the constraint force divided 
by the corresponding mass: 

p 

ZF c /m minimized. (Gauss' Principle) 

This principle provides a unique prescription for carrying out simulations of 
nonequilibrium flows. Several consequences of this principle have been explored. 

p 

Constant-energy and constant-temperature (with temperature T - p /3mk, where k 
is Boltzmann's constant) shear flows yield nearly identical viscosity coefficients 

Q 

even far from equilibrium. Self-diffusion, shear and bui k viscosity, and heat 
conduction can be studied in very small Systems (two particles for the first three 
transport coefficients, and three for heat flow). 10-12 When periodic bound- 
aries are used, to eliminate surface effects, the resulting two-particle transport 
coefficients lie within a factor of two of the many-body thermodynamic limit. 

Thus it appears that nonequi1 ibrium results have an intrinsic number-dependence no 
greater than the equilibrium dependence, 1/N, where N i s the number of particles. 1 

Gauss' principle can lead to erroneous results. If a simulation is carried 
out with fixed heat-flux vector and fixed temperature the corresponding thermal 
conductivity can be in error, in a dense fluid, by several percent. 14 

Nevertheless, a variety of simulations led to accurate values for the trans¬ 
port coefficients of simple fluids and solids. These new data suggested useful 
corresponding States relations between the transport coefficients and equilibrium 















These same calculations also suggest that the nonequilibrium conditions pre- 
vaillng in a strong shockwave be simulated by dynamical adiabatic compressions, 
with the time of the compression chosen to correspond to the picosecond timescale 
of rising stress in a real shockwave. In such a rapid compression it can be an- 
ticipated that coupling of the vibrational and transl ational modes is incomplete 
and that highly nonequilibrium energy distributions resulting from the lack of 
coupling produce Chemical reactions at far different rates than would be predicted 
from Arhenius kinetics. The relaxation process is just beginning to be studied in 
a serious way by using nonequi librium molecular dynamics.^ These results 
indicate that Nose's Dynamics is best suited to simulations far from equilibrium. 

B. Nose's Dynamics 

Nose pointed out that if the friction coefficient z in the generalized 
equation of motion 


mr - F(r) - z p , 

relaxes according to a first-order differential eguation 

z - a(p 2 -3mkT) , 

where a is fixed, the equilibrium distribution is identical to Gibbs' canonical 
distribution for the temperature T. 

This approach has been generalized to nonequi1ibrium Systems and appears to 
be very successful. Just as is the case with Gauss' principle it is possible to 
analyze Systems close to equilibrium and to show that Nose's dynamics provides 
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exact linear transport coefficients. Bymaking temperature depend upon space and 
time, and by introducing also stress as an independent (tensor) variable it is 
straightforward to generalize Nose''s approach to the treatment of nonequi 1 ibrium 
problems such as the adiabatic compression of hexanitrobenzene. These simulations 
have not yet been carried out, because they are relatively complex, but the tools 
for doing this are now available. The model for crystalline hexanitrobenzene i s 


described i n the next section. 









III. HEXANITR08EN2E NE 

Hexanitrobenzene, is a crystalline high explosive, about twice 

the density of water, first synthesized and characterized in Russia^. 

Initial attempts to make this material in the United States failed. Later, when 
it became clear that hexanitrobenzene is sensitive both to water and to light, 
military interest in it subsided. Nevertheless, from the scientific standpoint 
this substance is extremely interesting. Because it contains no hydrogen it is 
reasonable to use a classical description for its motion. 

The molecules have sixfold rotational symmetry, but wi th the nitro groups 
aligned like blades of a propellor to alleviate steric repulsion of neighboring 
oxygens. The crystal structure was described in 1966 by Akopyan, Struchkov, and 
Oashevskii. For reasons that are not clear the sixfold symmetry i s broken i n the 
crystalline phase. The molecules are arranged in layers with nearly hexagonal 
symmetry, with each nearly-equi1ateral triangle of molecules having a "long" side 
about one percent longer than its two "short" sides. These nearly hexagonal 
sheets are stacked, i n a nearly-face-centered-cubic arrangement: ABCABCABC..., 
with a spacing between planes 4.08A, considerably smaller than the lengths of the 
triangle legs, 9.05 and 9.13A. 

Discussions with quantum chemists indicated that a reliable simulation of the 
breakup of a hexanitrobenzene molecule i s beyond the State of the art. A likely 
mechanism i s the breakup, due to bending, of one of the six nitro groups. Faced 
with a lack of detailed information on the breakup mechanism we modelled the 
molecule as six carbon atoms and six nitro groups, with 6x3+6x3=36 degrees of 
freedom per molecule. Interactions within the molecules can be treated in the 
quasiharmonic approximation. Interactions between molecules were based on a 
truncated Lennard-Jones potential to which a Hooke's-law potential was added i n 
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such a way as to bring both the potential and its first derivative (minus the 
force) to zero at the potential cutoff. 

Exploratory calculations, in two dimensions, established that the energy was 
lowered by about 7 % i f the molecules were allowed to rotate. In these calcula- 
tions the Lennard-Jones potential was cut off at twice the distance to the poten¬ 
tial minimum which was in turn three times the C-C or C-NC^ separation distance. 

In the three-dimensional case both the nearly-hexagonal and nearly-face-centered 
forms of the structure were studied. The least-energy rotation, just over 8 
degrees, was approximately equal to that found i n two dimensions. The intermolec- 
ular and interplanar spacings were 8.75$ and 2.67$, respectively. The correspond- 
ing experimental values are (9.05, 9.05, 9.13)$ and 4.08$. The eiergies of the 
two forms studied agreed to four-figure accuracy. Because the range of the 
potential exceeds twice the interplanar spacing the numbers are in fact slightly 
different. 

A final investigation was designed to match the experimental spacings of 
hexanitrobenzene, 9.08$ and 4.08A. A stress-free unit cell was constructed using 

O O 

a potential minimum at 4.56A and a C-C + C-NO 2 separation of 2.79A. Again the 
two forms of crystal agreed i n energy to four-figure accuracy. A well depth of 
9.6xl0 -14 ergs was chosen approximately to reproduce the estimated heat of 
vaporization of the crystal. A detailed report of these calculations is in 

preparation'8. 





IV. RECOMMENOATIONS 

It i s recommended that a study of the normal-mode energi es, bond displace- 
ments, and angular distortions in the model of hexanitrobenzene be studied as a 
function of initial temperature, shockwidth and strength of shockwave. These 
calculations would be useful in assessing the validity of competing theories of 
vibrational relaxation, as described by Holian^ and in further developing 
technigues for predicting the sensitivity of new explosive Chemicals to detonation 
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